How much does it matter and can we understand it with MCMC
Using wikipedia and others names below: https://en.wikipedia.org/wiki/Dead_time
Bottom work is based on Adams RP, Murray I, MacKay DJC. Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities. Proceedings of the 26th Annual International Conference on Machine Learning; Montreal, Quebec, Canada. 1553376: ACM; 2009. p. 9-16.
In [1]:
%matplotlib inline
from pprint import pprint
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as mc
import spacepy.toolbox as tb
import spacepy.plot as spp
import tqdm
from scipy import stats
import seaborn as sns
sns.set()
matplotlib.pyplot.rc('figure', figsize=(10,10))
matplotlib.pyplot.rc('lines', lw=3)
matplotlib.pyplot.rc('font', size=20)
%matplotlib inline
In [2]:
np.random.seed(8675309)
Rate1 = 10
Rate2 = 20
deadtime = 0.05
times1 = np.random.exponential(1/Rate1, size=1000)
times2 = np.random.exponential(1/Rate2, size=1000)
times = pd.DataFrame({'between1':times1, 'between2':times2})
times['times1'] = np.cumsum(times['between1'])
times['times2'] = np.cumsum(times['between2'])
# get the times of bunch of hits
In [3]:
fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(20,8), squeeze=False)
print(ax[0][0])
times[['between1', 'between2']].plot(ax=ax[0][0])
ax[0][0].set_ylabel('Time between')
times[['times1', 'times2']].plot(ax=ax[0][1])
ax[0][1].set_ylabel('Wall time')
plt.tight_layout()
plt.figure(figsize=(20,8))
times[['between1', 'between2']].hist(bins=30, figsize=(20,8))
Out[3]:
In [4]:
# make cts per unit
realcps1 = np.histogram(times['times1'], bins=np.arange(0, times['times1'].max()))[0]
realcps2 = np.histogram(times['times2'], bins=np.arange(0, times['times2'].max()))[0]
In [5]:
# then cull out all points and are closer than deadtime to each other
times['between1dt'] = times['between1'].copy()
times['between1dt'].loc[times['between1dt'] < deadtime] = np.nan
times['between2dt'] = times['between2'].copy()
times['between2dt'].loc[times['between2dt'] < deadtime] = np.nan
times['times1dt'] = np.cumsum(times['between1dt'])
times['times2dt'] = np.cumsum(times['between2dt'])
In [6]:
times[['times1', 'times1dt']].plot(figsize=(10,7))
times[['times2', 'times2dt']].plot(figsize=(10,7))
Out[6]:
In [7]:
dtcps1 = np.histogram(times['times1dt'].dropna(), bins=np.arange(0, times['times1dt'].dropna().max()))[0]
dtcps2 = np.histogram(times['times2dt'].dropna(), bins=np.arange(0, times['times2dt'].dropna().max()))[0]
In [8]:
df = pd.DataFrame({'dtcps1':dtcps1, 'realcps1':realcps1[:len(dtcps1)]})
df.plot()
# TODO looks wrong!!
Out[8]:
In [ ]:
In [9]:
np.random.seed(8675309)
size = 100
n_rate = 6
data = np.empty((n_rate, size))
for ii, d in enumerate(tb.linspace(1, 100, n_rate)):
data[ii] = np.random.poisson(d, size=size)
print(data.shape)
In [10]:
plt.plot(data.T);
From the reference:
The relationship of the real counting rates with the measured counting rates is well known for these two basic models. If we denote by N the real counting rate, by M the measured counting rate, and by τ the dead time and considering that the non-dead time disturbed distribution is Poissonian, M and N are related by
$N=\frac{M}{1-M\tau}$
for the case of a non-paralyzable dead time, and by
$M=Ne^{-N\tau}$
In [11]:
with mc.Model() as model0:
mu = mc.Uniform('mu', 0, 1e6, shape=n_rate)
dat = mc.Poisson('dat', mu=mu, observed=data.T, shape=n_rate)
# det2 = mc.Poisson('d2', mu=mu, observed=d2[0:10])
start = mc.find_MAP()
trace0 = mc.sample(10000, start=start)
In [12]:
mc.summary(trace0)
In [13]:
ax = mc.traceplot(trace0, lines={'mu':tb.linspace(1, 100, n_rate)})
In [ ]: